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ABSTRACT 

We follow numerically the nonlinear evolution of the Parker instability in the presence of 
phase transitions from a warm to a cold Hi interstellar medium in two spatial dimensions. 
The nonlinear evolution of the system favors modes that allow the magnetic field lines to 
cross the galactic plane. Cold Hi clouds form with typical masses ~ 10 5 M , mean densities 
~ 20 cm -3 , mean magnetic field strengths ~ 4.3 /iG (rms field strengths ~ 6.4 /iG), mass-to- 
fiux ratios ~ 0.1 — 0.3 relative to critical, temperatures ~ 50 K, (two-dimensional) turbulent 
velocity dispersions ~ 1.6 km s , and separations ~ 500 pc, in agreement with observations. 
The maximum density and magnetic field strength are ~ 10 3 cm -3 and ~ 20 /iG, respectively. 
Approximately 60% of all Hi mass is in the warm neutral medium. The cold neutral medium 
is arranged into sheet-like structures both perpendicular and parallel to the galactic plane, 
but it is also found almost everywhere in the galactic plane, with the density being highest 
in valleys of the magnetic field lines. 'Cloudlets' also form whose physical properties are in 
quantitative agreement with those observed for such objects by Heiles (1967). The nonlinear 
phase of the evolution takes < 30 Myr, so that, if the instability is triggered by a nonlinear 
perturbation such as a spiral density shock wave, interstellar clouds can form within a time 
suggested by observations. 
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1 INTRODUCTION 

Understanding the formation of interstellar clouds and cloud com- 
plexes is a problem of longstanding importance, one which has re- 
ceived increased attention in recent years. Since Hi clouds are the 
progenitors of molecular clouds, the stellar birthplaces, the process 
by which they form is of fundamental importance to galactic evo- 
lution as well as to the early stages of star formation. 

Proposed mechanisms for interstellar cloud formation 
fall into three categ or ies: stochastic c oagula ti on of s maller 
clouds JOortl 119541; iField & Saslawl 1 19651 : iKwanl 1 19791 : 
IScoville & Hershl I1979D. colliding streams of turb ulent flows 
jBallesteros-Paredes. Hartmann & Vazquez-Semadeni| 19991: 



Hartmann. Ballesteros-Paredes & Berginl 1200 ll : iHeitsch et al.l 
2006 ) , or collective effects involving instab i lities dMouschoviasI 
19741 : iMouschovias, Shu & Woodward! 1 1974 iBlitz & Shulll98(ifh 
Stochastic coagulation is not feasible because collisional ag- 
glomeration proceeds too slowly and, also, collisions are more 
likely to result in disruption rather than merger. Turbulence-driven 
formation of atomic-hydrogen clouds is problematic because no 
realistic mechanism has yet been demonstrated to sustain the 
required large-scale, ordered, focused motions for ~ 10 Myr. 
In addition, the short lifetimes and highly supercritical mass- 
to-flux ratios of molecular clouds thought to form within these 
turbulent 'atomic precursors' are in conflict with observations 



dMouschovias. Tassis & Kunzl l2006h . Q By contrast, large-scale 
dynamical instabilities can amass material over lengthscales much 
greater than the sizes of diffuse atomic-hydrogen clouds in a time 
consistent with observa tional constraints. 

In the iParkea i 1966ft instability, undulations of initially hori- 
zontal (i.e., parallel to the galactic plane) magnetic field lines allow 
matter to slide down along field lines under the action of the ver- 
tical component of the galactic gravitational field. The raised por- 
tions of the field lines are thereby relieved of the confining weight 
of the gas and rise even higher, until the tension of the field lines 
increases sufficiently for a final equilibrium state to be established 
( Mous chovias!! 19741) . Mouschovias et al. (1974) have argued that 
the Parker instability is triggered in spiral arms behind spiral den- 
sity shock waves, leading to the formation of cloud complexes in 
the valleys of the curved magnetic field lines. 

Several consequences of this mechanism for the formation of 



The typical separation between the first appearance of CO and the 
peak of Ha in spiral arms corresponds to a timescale ~ 10 7 yr. 
This effect cannot possibly be due to dus t obscuration, as propos ed in 
i Garcia-Burillo. Guelin & Cernicharol < 19931) . |Kim & Ostrikeil <2006l) . and 
lElmegreenl 420071) . A similar shift occurs between the first appearance of 
CO and the peak o f Pact, which does not suffer from dust obscuration 
IScoville et alj200ll) . 
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interstellar clouds are in excellent, in fact unique, agreement with 
observations: 

(i) The implosion of individual clouds by shock waves within 
these complexes leads to the formation of OB associations and gi- 
ant Hll regions aligned along spiral arms 'like beads on a string' 
and separated by regular intervals of 500 — 1000 pc, as observed 
both in our galaxy and in external galaxies jWesterhoutlll963l ; lKerJ 
1 19631: iHodgell 19691 : lMorgarJll970h . 

(ii) The observed moderate enhancement in synchrotron 
radiation from the interarm to the arm region of spiral 
galaxie s l Mathewson. van der Kruit & Brouwll 19721 ; I van der Kruit 
1973a.hc) found a natural explanation in the context of the Parker 
instability (Mouschovias et al. 1974; Mouschovias 1975a). The re- 
lief of magnetic stresses by buckling of the field lines in the ver- 
tical direction in spiral arms leads to an increase in the magnetic 
field strength that is much smaller than that expected from one- 
dimensional compression. Moreover, the tendency for the cosmic- 
ray pressure to distribute itself uniformly along magnetic flux tubes 
means that, after the instability develops, most of the relativistic 
particles are found in the regions where the field is weak, i.e., in 
the inflated portions of the flux tubes. In other words, the cosmic- 
ray and magnetic-energy densities are not in equipartition on scales 
< 1 kpc; there is an anticorrelation between them. 

(iii) The anticorrelation between the cosmic-ray density on the 
one hand and the magnetic-field strength and gas density on the 
other hand also implies that the intensity of synchrotron radiation 
along spiral arms should not vary significantly, as observed. This 
is in sharp contrast with expectations based on correlation between 
the field strength and the cosmic-ray density (e.g., equipartition be- 
tween magnetic and cosmic-ray energy densities). 



Numerical simulations of the Parker instabil- 
ity with applications to the interstellar medium have 
been 



carried out by iBasu. Mouschovias & P aleologou 
1997: hereafter BMP97I), iKim et alJ dl998h , iKim. Ryu & Jonesl 



20011) , and lKim, Ostriker & Stonej (2002). Several questions con- 



cerning the nonlinear evolution of the instability were addressed. 
For example, which mode dominates the nonlinear growth of the 
instability on timescales of interest in the interstellar medium 
(ISM); how fast the instability grows in the nonlinear regime; and, 
how great a density enhancement is generated in the galactic plane 
by downflow along the magnetic arches. A common finding in 
these papers is that the Parker instability by itself cannot be the 
agent responsible for interstellar cloud formation; the maximum 
density enhancement during the nonlinear evolution is only a 
factor of ~ 2. Instead, the Parker instability may be thought of as 
setting the stage upon which smaller-scale processes, such as cloud 
formation, star formation, and subsequent supernova explosions, 
act out their individual roles and collectively contribute to the 
structure and evolution of the ISM. However, Mouschovias (1974) 
and BMP97 also noted that the Parker instability may be directly 
responsible for the formation of GMCs, or cloud complexes, if 
it can generate a sufficient den sity enhancement near the galactic 
plane for a critica l press ure (Field, Goldsmith & Habin 4 1 19691 : 
iKapl an & P ikel' ne 3 1 19701) to be exceeded, thereby initiating a 
phase transition to a cooler, denser cloud ph ase. This point has 
been re-emphasized bv lParker & Jokipii (200(3). 

We follow numerically the nonlinear evolution of the Parker 
instability in two spatial dimensions accounting for phase transi- 
tions. In Section [2] we present the basic equations and physics of 
the instability. In Section [3] the dimensionless equations are given 
and the numerical method of solution is described. The results of 



the nonlinear development of the instability with phase transitions 
from a warm neutral medium (WNM) to a cold neutral medium 
(CNM) are presented in Section|4] followed by a summary and dis- 
cussion in Section [5] Contact between theory and observations is 
emphasized. 



2 FORMULATION OF THE PROBLEM 
2.1 Basic Equations 

We consider a conducting fluid of mass density p, thermal pres- 
sure P, and temperature T, threaded by a frozen-in magnetic field 
B, in an external gravitational field g. The magnetohydrodynamic 
(MHD) equations describing this system are 



| + v. M = o, 



9(pv) 
dt 



+ V ■ (pvv) 



(la) 



-VP + pg + — (V X B) x B, (lb) 

47T 



(lc) 



§ =VX(„XB), 

fin 

^ + V • (uv) = -PV • v - pL + V • (kVT) , (Id) 

where v is the velocity and u is the internal energy density. The 
thermal pressure is a function of both density and temperature, 
P = P(p, T), and is determined by a balance of heating and cool- 
ing mechanisms in the ISM. It is obtained from the internal energy 
density of the gas through the ideal gas relation, P = (7 — l)it, 
where 7 is the ratio of specific heats (= 5/3). In equation j kit , 
C = pA — r is the loss function (i.e., net cooling rate per unit 
mass), where A = A(p, T) and T are the cooling and heating rates, 
respectively; n is the thermal conductivity (see eq.U2lbelow). 

We use an equation of state that allows the existence of two 
phases, a warm one with T ~ 10 4 K and density n ~ 0.1 — 1 cm" 3 



and a cold one with T ~ 10 K and n ~ 10 - 100 cm 



Follow- 



ing [Nagashima^Jnutsuka^^byam^ |2006), we adopt the follow- 
ing heaiing^idcooling^ates0 



firriH- 



pmyT 
A(T) 



2 x 10" 26 erg s" 1 . 



= 10 exp 



-1.184 x 10^ 
T + 1000 



+ 1.4 x 10~ 2 \/Texp ( ) cm 3 



(2a) 



(2b) 



where p — 1.27 is the mean mass per particle in units of the 
atomic-hydrogen mass, accounting for 10% He abundance by num- 
ber. The minimum and maximum pressures for the coexistence of 
the warm/cold phases in a static equilibrium are P mm /fcB — 1600 
Kcm -3 and Pmax/^B — 5000 K cm" 3 . The corresponding transi- 
tion temperatures T m i n = 185 K and T max = 5012 K define cold 
(T < T m i n ), warm (T > T max ), and intermediate temperature 
(Tmin < T < Tmax) phases. Since we do not include the addi- 
tional phase transition to molecular hydrogen in our loss function, 
we enforce isothermality for any gas element whose number den- 
sity exceeds rii so = 60 cm" 3 . The corresponding temperature of 
this isothermal branch is Ti so = 51.3 K. The thermal equilibrium 
curve implied by these heating and cooling functions is shown in 



2 The function al form of these heating and cooling rates has its origin in 
earlier work bv lPenston & Brownl fT970) and references therein. 
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Figure 1. Thermal equilibrium curve (solid line) implied by the heating 
and cooling functions given by equations [2a| and [2b] The dashed lines are 
isotherms marking the transition temperatures between different phases of 
the ISM; the WNM corresponds to T > T max = 5012 K and the CNM 
to T < T m ; n = 185 K. The units of P and p are their initial values 
(in the 'reference', stratified state) at the galactic plane 2 = 0. Typically, 
Prcf(0) = 2.1 X 10 -24 g cm -3 , corresponding to a number density 
= 1 cm -3 , and C rc f(0) = 5.73 km s _1 , corresponding to a temperature 
T rcf (0) = 5012 K (see § 3.1). 

Figure[T] The pressure (ordinate) and density (abscissa) are normal- 
ized to their initial midplane values (see § |3.1| below). The dashed 
lines denote the T max -, T m i n -, and Ti so -isotherms. 

Note that an investigation of the Parker instability account- 
ing for phase transitions is not the same as an investigation of the 
therm al instability occurring alongside the Parker instability dAmesI 
Il973h . The lengthscales involved in the two instabilities are so rad- 
ically different (~ 1 pc for the thermal instability and ~ 1 kpc 
for the Parker instability) that the two instabilities proceed in paral- 
lel without much effect of one on the other. The key element in the 
possibility of the Parker instability leading directly to the formation 
of massive Hi clouds and Hi cloud complexes is not so much the 
detailed nature of the thermal instability but the fact that a tran- 
sition from the warm to the cold phase can take place by some 
means (e.g., even if the equation of state on a log P - log p dia- 
gram had only a horizontal branch joining the two phases, instead 
of the branch with a negative slope). 



3 METHOD OF SOLUTION 

3.1 Initial and Boundary Conditions 

We consider a local Cartesian coordinate system (x,y,z) corre- 
sponding to cylindrical polar coordinates (r, if), z) in a galactic disk, 
with the range of the x (radial) coordinate being <C R (the distance 
from the galactic center) so that curvature terms can be ignored. 
The magnetic field is initially in the azimuthal direction (y), which 
we take to be along a spiral arm. Because the motions involved in 
the Parker instability are predominantly azimuthal, rotation can be 
ignored (see Shu 1974; also Mouschovias et al. 1974). 

Observations show that the interstellar magnetic field is pre- 
dominantly parallel to the galactic plane, although arches of mag- 
netic field lines rising hi g h above the plane are also revealed 
dMathewson & Fordlll97rj|) . |ParkeJ Jl966h assumed a zeroth-order 
initial (reference) equilibrium state, which, for simplicity, has field 
lines exactly parallel to the galactic plane (z — 0), so that B ro f = 
B le f(z)e y . The vertical galactic gravitational field (due to stars) is 



taken to be constant and to reverse its direction across the galac- 
tic plane; i.e., g = — g(z)e z , where g(z) = — g{— z) = g = 
constant > 0. More comp licated gr a vitatio n al fields ha v e bee n 
considered analyti c ally b y iKellmanl fl972h, iGiz & Shul dl993h. 



iKim. Hong & Rvul dl997h. iKim & Hond Jl998h . and 



Kim et al 

(2000) and numerically by Santillan et al. (2000) and Fran co et al 



120021) . The main effect of a non-constant gravitational accelera- 
tion is to increase the growth rate of the instability and shorten the 
separation of the magnetic valleys. 

Although we regard the Parker instability as a driven (or 
forced) instability in galactic spiral arms and, therefore, the 'ini- 
tial' state is not an equilibrium one, we nevertheless use Parker's 
initial equilibrium state as a reference state, from which we obtain 
the input parameters of the problem. This approach ensures that (i) 
any 'final' outcome of the evolution of the gas-field-gravity system 
can be reached from Parker's stratified equilibrium state through 
continuous deformation of the field lines; and, (ii) the results of 
the present calculation can be compared and contrasted with those 
of past calculations on the nonlinear evolution of the Parker insta- 
bility in the ISM. We thus take the ratio of magnetic and thermal 
pressures, 



a = B Icf /8nP le{ . 



(3) 



to be constant in the reference equilibrium state, in which the den- 
sity is given by 



p le f(z) = p re f (0) exp 



where 



H=(l + a)Cl { (0)/g 



M 

H 



(4) 



(5) 



is the total scale height of the gas in the reference state and C re f (0) 
is the isothermal sound speed in the reference state at the galactic 
plane. With T rcf (0) = 5012 K, g = 3 x 10" 9 cm s -2 , and a ~ 1, 
a value consistent with observations, one finds that H ~ 71 pc. 
The magnetic field strength in this reference state is given by 



B ref (z) =S lef (0)exp( 



(6) 



The corresponding pressure and temperature distributions are de- 
termined by requiring the internal energy to satisfy local ther- 
mal equilibrium, based on the heating and cooling functions 
given in Section 12.11 For simplicity, we ignore the effect of 
cosmic rays in this study, although their effects on the evo- 
lution of the Par ker instability are w e ll known by both an- 



alytic arguments dMouschoviasI [l975bl. 1 19961; lHanasz & Leschl 

33j;L 



J and numerical work faanasz & Leschl2003T ; Ryu et al .120031 ; 
iKuwabara. Nakamura & Koll2004l) . Mainly, they tend to speed up 
the instability. 

The solution of the linear instability problem can either have 
field lines with reflection symmetry about z = (midplane- 
symmetric or 'odd' mode) or allow field lines to cross z — 
(midplane-crossing or 'even' mode). In the latter case, the quan- 
tities p, B, and v are invariant with respect to a reflection about 
z = followed by a translation in the y-direction by an amount 
±A H /2, where X y is the horizontal (along the field lines) wave- 
length of the perturbation. (The terms 'even' and 'odd' refer to the 
algebraic form of the velocity perturbation imposed on the refer- 
ence equilibrium state.) The 'odd' mode in an isothermal system 
does not lead to any appreciable density enhancement, and is there- 
fore less likely than the 'even' mode to produce interstellar clouds 
in a nonisothermal system. Moreover, BMP97 found that the 'even' 
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mode is naturally selected by the nonlinear evolution of the system. 
Hence, we study a model with an initial perturbation that has, and 
at later times leads to evolution that preserves, the symmetry of the 
even mode[f]This symmetry allows us to study a region of horizon- 
tal extent Y = A a /2. 

The equations are put in dimensionless form by choosing 
units natural to the problem. The units of velocity [v], density 
[p] , and magnetic field strength [B] are, respectively, the values of 
the isothermal sound speed C rc f (0), density p rc f (0), and magnetic 
field B TC { (0) at the galactic plane|j The unit of length [L] is the 
thermal-pressure scale height Cf cf (0)/g. The implied unit of time 
[t] is C ro f (0)/g, the sound-crossing time across a thermal-pressure 
scale height. While the evolution of the Parker instability depends 
on only one free parameter, a, specific values of p re f (0), g, and 
T re f (0) must be chosen to put the loss function C in dimensionless 
form. We choose p rc f(0) = 2.1 x 10 -24 g cm" 3 (correspond- 
ing to a number density of 1 cm -3 ), g — 3 x 10~ 9 cm s -2 , and 
T re f (0) = 5012 K. (This corresponds to a thermally stable state in 
the warm phase.) Choosing the free parameter a = 1, we find the 
following values for the units: [v] = 5.73 km s~ , [L] = 35.45 pc, 
[t] = 6.05 Myr, and [B] =4.17 pG. 

We choose Y = 12 X [L] = 425.4 pc, so that X y ~ A^max 
(the horizontal wavelength corresponding to the maximum growth 
rate of the linear instability). The density and velocity are assumed 
to have reflection symmetry about the planes y = and y = Y, 
where the field lines remain locally horizontal. The upper and lower 
boundaries, where the field lines remain undeformed, act as 'lids' to 
the system, and are placed at z = ±Z — ±25 x [L] = ±886.25 pc, 
so that they are far enough from the galactic plane not to affect the 
evolution of most of the matter between them0 (In the reference 
state the density at z = ±Z decreases to exp( — 12.5) = 3.73 x 
10 -6 of its value at the galactic plane.) The velocity perturbation is 
taken to have the form 



Sv z (y,z) 



-eCr 



(5) 



(7) 



where e = 0.2, so that the kinetic energy of the perturbation is 
less than 1% of any other form of energy (gravitational, thermal, 
or magnetic) in the system. This choice of the amplitude of the 
perturbation does not, of course, represent the highly nonlinear dis- 
turbance introduced in nature by a spiral density shock wave, and it 
also has the consequence that the instability takes a longer time to 
enter the nonlinear regime than it does in nature. This choice, nev- 
ertheless, allows us to follow the development of the linear phase 
of the instability as well as avoid the introduction of an arbitrari- 
ness in the problem by using such large perturbations that dominate 
from the outset the system's evolution. To mimic the natural phe- 
nomenon of the instability driven by a spiral density shock wave 
without having to resort to a three-dimensional calculation, we set 



3 The simulations by Kim et al. (1998, 2001) could only study the 'odd' 
mode, due to their assumed symmetry about the galactic plane. BMP97, 



Santill an et al U2000|) . and Kim et al. (2002) considered both the 'even' and 
'odd' modes. 

4 One could choose the unit of B to be [o87rp ref (0)C 2 cf (0)] 1 / 2 , but it is 
more convenient to use B IC { (0); the two are identical for a = 1. 

5 Several numerical studies of the Parker instability miss t his crucial point. 
For ex ample, the simulations by Kim et al. (2002) and iKim & Ostrikej 
120061) have a total vertical box size of only 8H (with a = 1/2). Our 
total vertical box size, by contrast, is ~ 25-ff ~ 1.8 kpc. A small verti- 
cal size of the comp utational box suppresses the ins tability [see eq. (32) of 

r by [Mouschovias i ll99d> l. 



MouschoviaJ ll974l) ; also, eq. (22) of the 1 



the origin of time (t = 0) at the instant the linear phase of the evo- 
lution has led to an increase of the maximum density in the system 
by a factor of 2. Times earlier than this instant are assigned negative 
values in what follows. 



3.2 The Dimensionless Equations 

We denote the dimensionless time by r and the dimensionless co- 
ordinates (y, z) by (y' , 2'), so that 

d g d 



and 



9 



V. 



8t C ref (0) dl 

All other dimensionless quantities are adorned with a tilde. The 
dimensionless equations are then 



dp 



+ V'-(pv) = 0, 



(8a) 



^^±V' ■ (fivv) = -V'P-pe z -aV'B 2 +2aB ■ V'B , (8b) 

^ = V'x(JxB), 
or 

da 



— + V' • (uv) = -PV' -v - pL + V' • (kV'T) , 



(8c) 
(8d) 



where P = (7 — l)u = pT is the dimensionless pressure, and 

u ~ C _ _ K<?T ro f(0) 



C 



p ref (0)C7 2 cf (0) ' 5 C7 re f(0) ' Prcf (0)C r 5 of (0) 



, (9) 



are the dimensionless internal energy, loss function, and conduc- 
tivity, respectively. The dimensionless density, magnetic field, and 
velocity perturbations in the reference state are, respectively, 



Prcf (2') = exp ( - 



B ie f(z') = exp 



1 + al ' 



2(1 + a) J ' 



and 



~ , , ,, , \ nz 

MV ,2 ) = -e cos ( — 1 cos I — 



(10a) 
(10b) 

(10c) 



These equations are subject to the boundary conditions 

v z (y',±Z') = 0, v y {0,z')=0, v y (Y',z')=0, (11a) 
B z (y',±Z')=0, B z (0,z')=0, B z {Y',z')=Q, (lib) 



dp(y',z" 



dy> 
dv z {y',z') 



= 0, 



y'=0,Y' 



dy> 
dB y (y',z') 



dz' 
dv y (y',z' 



y'=o,Y< 



8y' 



= 0, 



dz' 
dB y (y',z') 



= 0, (11c) 



0, (lid) 



V'=0,Y> 



dz 



= 0. (lie) 



>=±z> 



3.3 Numerical Issues 

We integrate equations ( f8ab - J8db using a version of the publicly- 
available Zeus-MP code jHaves et al.l2006t) distributed by the Lab- 
oratory for Computational Astrophysics of the University of Cali- 
fornia at San Diego. Zeus-MP uses a time-explicit, operator-split, 
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Figure 2. Time evolution of the (dimensionless) central (a) density p c , (b) temperature T c , and (c) magnetic field strength B c at (y' , z') = (0, 0), the bottom 
of a magnetic valley. The units of p c , T c , and B c are their values at the galactic plane, z' = 0, in the initial reference state, namely, n re f (0) = 1.0 cm -3 , 
^rcf (0) = 5012 K, and B rc { (0) = 4.17 ^G. The unit of time is the sound-crossing time across a thermal-pressure scale height in the reference state, i.e., 
C re f(0)/ a =6.05 Myr. 



finite-difference method for solving the MHD equations on a stag- 
gered mesh, capturing shocks via quadratic and linear artificial vis- 
cosities. De t ails c oncerning the algorithms used can be found in 
lHaves et al.l bOOfjh and references therein. 

We have made modifications to this code in order to follow 
phase transitions. Because cooling times can be very short com- 
pared to the timescales relevant to cloud formation, the energy 
equation update from the net cooling terms is solved implicitly us- 
ing Newton-Raphson iteration. At the start of each iteration, the 
timestep is initially computed from the Courant-Friedrichs-Levy 
condition using the sound speed, Alfven speed, and conduction 
parameter (see below). This is followed by a call to the cooling 
subroutine. The change in temperature within each zone is lim- 
ited to 10% of its initial value. If this requirement is not met for 
all cells in the computational domain, the timestep is reduced by 
a factor of 2, and the imp licit energy update is recalculated. As in 
IPiontek & Ostrikerl d2004) . we have found that the timestep is re- 
duced at most once or twice per convergent timestep. 

The update from the conduction operator is solved explic- 
itly, using a second-order a ccurate differencing of V • (kVT). As 
iKovama & Inutsuka ((2004) pointed out, the importance of incorpo- 
rating conduction in simulations that contain thermally unstable gas 
has sometimes been overlooked in past work. Without conduction, 
the growth rates for thermal instability are greatest at the smallest 
lengthscales, and unresolved growth at the grid scale may occur. 
The inclusion of conduction has a stabilizing effect on the thermal 
instability at small scales, and the conduction parameter can be ad- 
justed to allow spatial resolution of the thermal instability on the 
computational grid. We tune the thermal conduction coefficient k 
so that the number of zones in a Field length satisfies Af / Ay = 4 
in the thermally unstable regime|f] 



i - 



dlnA 
dlnT 







(12) 



As a parcel of gas transits (thermodynamically) from an unstable 
state to a stable state, the thermal conductivity vanishes in a contin- 
uous fashion. 

The timestep is necessarily limited not only by the usual sound 
and Alfven cell-crossing times, but also by a diffusion timescale 
that is set by the thermal conduction parameter: 



(At) c 



(Ay) 2 



4( 7 - 



(13) 



The quadratic dependence on the grid spacing is a general feature 
of any numerical problem that includes diffusion. This dependence 
is particularly restrictive when high resolution is required to accu- 
rately capture important physical processes. In this work, a high 
resolution is needed so that the thermal conduction parameter is 
restricted to reasonable values. 

The problem is solved on an orthogonal, nonuniform grid. The 
horizontal (y) grid is spatially uniform, with 2048 zones and a res- 
olution Ay ~ 0.21 pc (i.e., Ay' ~ 5.9 x 10~ 3 ). The vertical (z) 
grid contains 4096 zones and is in general nonuniform, with the 
resolution being greatest near the galactic plane. There are 2048 
uniformly-spaced zones in the 15 initial thermal-pressure scale 
heights straddling the galactic plane, so that Az ~ 0.26 pc (i.e., 
Az' ~ 7.3 x 10" 3 ); the remaining 2048 zones are logarithmically- 
spaced above and below this uniform-grid region. The transition 
between the uniform and nonuniform regions is chosen so that the 
grid spacing varies smoothly. 



6 The Field length is A F = 2tt{(p 2 A/kT)[1 - (din A/dlnT)]}- 1 / 2 
when A is a function of temperature. It is essentially the critical lengthscale 
inside which conduction stabilizes the thermal instability. 



4 RESULTS 

In Figure[2] we show the time evolution of the central (dimension- 
less) (a) density, (b) temperature, and (c) magnetic field strength 
at (y' , z') — (0, 0), which is located at the bottom of a magnetic 
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T = -13 T = -8 T = -3 




Figure 3. Contour plots of the density (color) and magnetic field lines (solid lines) at times r = —13, —8, —3, 2, 7, and 12 in the (y' , z' ) plane. The colorbar 
gives the correlation between the colors and the natural logarithm of the density [e.g., '+4' denotes a dimensionless density of exp(+4) ~ 54.6]. The unit of 
length is the thermal-pressure scale height in the reference state, (Q)/g = 35.45 pc, the unit of density is p re f (0) = 2.1 X 10 — 24 g cm -3 , and the unit 
of time is C* rcf (0)/g = 6.05 Myr. 
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valley. As found by previous studies of the isothermal Parker insta- 
bility, there are three distinct phases of evolution: linear, nonlinear, 
and relaxation into a 'final' equilibrium state. 

During the linear phase (— 13 ^ r ^ 0), the velocities remain 
subsonic and the central density, temperature, and magnetic-field 
strength vary only slightly. Soon after the onset of the nonlinear 
phase (at r = 0), the central density, temperature and magnetic- 
field strength change very rapidly. After the density has increased 
by a factor of ~ 2, the pressure of many fluid elements has also 
increased above the critical value for a phase transition to occur 
to a cold (cloud) phase. This leads to a rapid density enhancement 
by more than two orders of magnitude in a time At < 3 (i.e., 
At < 18 Myr). During the phase transition, approximately isobaric 
evolution causes the temperature to decrease by about two orders 
of magnitude. The magnetic field strength increases only slightly 
during the early part of the nonlinear phase, since gas motions are 
primarily along field lines. However, after enough gas has accu- 
mulated in the valleys of the field lines, its weight compresses the 
field lines and the maximum magnetic-field strength increases by a 
factor of three (from ~ 4 /j,G to ~ 12 fiG). During the subsequent 
relaxation phase, the compressed gas rebounds because it had over- 
shot its equilibrium state, and eventually settles in a state that has 
a central density ~ 200 cm -3 . The magnetic field strength contin- 
ues to increase as more mass elements enter the cloud phase and 
the galactic gravitational field (i) compresses the cloud perpendic- 
ular to the galactic plane and (ii) flattens it in the galactic plane. 
The resulting physical conditions are such that, if we had included 
the gas' self-gravity in the calculation, individual self-gravitating 
clouds would separate out. The density oscillations evident after 
r ~ 4 occur about a 'final' equilibrium state; they are the result of 
motions left over by the cloud formation process. 

As explained in Section 3.1, the linear phase of evolution is 
expected to be very short-lived in nature because the instability is 
externally driven by a spiral density shock wave (see Mouschovias 
et al. 1974). For this reason, the 'origin' of time is chosen to be the 
instant by which the central density has increased by a factor of 2. 
It is thus seen that the time required to form cold interstellar cloud 
complexes from the WNM does not exceed ~ 30 Myr. 

Figure[3]exhibits contour plots of the density (color) and mag- 
netic field lines (solid lines) at times r = — 13, —8, —3, 2, 7, and 
12. For ease in visualization, we show a region that covers one full 
horizontal wavelength, —Y' ^ y ^ +Y', of the instability. Re- 
gions of high (low) density are colored red (blue) with a color con- 
tinuum between them (as shown in the colorbar to the right of the 
figures). By the end of the linear phase of the evolution (r = 0), 
the magnetic field lines have begun to deform significantly and the 
gas has begun its rapid descent into the magnetic valleys. The maxi- 
mum velocity has exceeded the sound speed by a factor of 1 .6 1 . The 
evolution up to this point is nearly indistinguishable from the evo- 
lution of the isothermal Parker instability (determined by BMP97), 
because few gas elements have transitioned into the stable CNM 
phase by this time (see below). By r = 2, the system is well into the 
nonlinear stage of its evolution. The field lines have become con- 
siderably arched, even near z = 0. There is a strong flow of gas 
[maximum velocity 2.63C rc f (0)] into the magnetic valleys, lead- 
ing to the formation of shocks. By r = 7, the rising portions of 
the field lines have almost stopped their expansion and are nearly 
at their 'final' equilibrium positions. 

In Figure [4] we display the large-scale velocity (arrows) and 
density (solid lines) fields at time r = 12. There is considerable 
density structure in the galactic plane. The velocity field shows 
prominent and distinctive S-shaped shocks extending far above and 



y max = 2.74C ref (0) 




Figure 4. Contour plot of the density (solid lines) at time r = 12 in the 
(y' , z') plane. The velocity vectors (arrows) are also shown and the max- 
imum value of the velocity [= 2.74C ro f (0) ~ 15.7 km s^ 1 ] is printed 
above the figure. The unit of length is (7 2 f (0) /<) = 35.45 pc. Note the 
presence of distinctive S-shaped shocks extending far above and below the 
galactic plane. 

below the galactic plane. Shocks of this character have also been 
see n in the BMP97 sim ulations of the isothermal Parker instabil- 
ity. iMouschoviasI i ll 9961) noted that, if such shocks are observed, 
they cannot possibly be mistaken for supernova remnants and can 
therefore be taken as evidence for the nonlinear development of the 
Parker instability in the interstellar medium. 

Figure^ shows the innermost —1 ^ y' ^ +1, —1 ^ z' ^ 
+ 1 part of the system shown in Figure [5] at time r = 12. The 
CNM has a clear sheet-like morphology extending high above the 
galactic plane, but it also has a signific ant extent in the galactic 
pla ne, in agreement with ob servations bv lHeiles & Jenkins! d 19761) 
an d lHeiles & Trolandl J2003I) . The maximum velocity in this region 
is Umax = 1.49C rc f (0) — 8.5 km s _1 , also in agreement with ob- 
servations. Gas elements slide down the deformed magnetic field 
lines and undergo a phase transition, from the WNM to the CNM. 
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Figure 5. Contour plots of the density (color) and magnetic field lines (solid lines) at time r = 12 in the (y' , z') plane, zoomed in on (a, left) a CNM 'sheet' 
in a magnetic valley and (b, right) a CNM cloudlet. In (a) the velocity vectors (arrows) are also shown and the maximum velocity [= 1.49C re f (0) ~ 8.5 km 
s _1 ] is given. In (b) the maximum density [= 9.80p r cf (0) ~ 2 X 10 — 23 g cm -3 ~ 10 cm -3 ] of the cloudlet is given. The unit of length is C 2 c j (0)/g = 
35.45 pc. 




Figure 6. Contour plots of the density (color) and magnetic field lines (solid lines) in the (y' , z') plane, zoomed in on a CNM 'sheet' in the galactic plane, 
at times (a, left) r = 7 and (b, right) r = 12. The velocity vectors (arrows) are also shown and are normalized to the maximum speed in the region [= (a) 
0.98C rc f (0) ~ 5.6 km s _1 ; (b) 0.27C re f (0) ~ 1.5 km s -1 ]. The unit of length is C 2 ef (0)/g = 35.45 pc. Note that the z' axis is stretched out of proportion 
in order to show the vertical structure of the sheet. 



A shock front marks the CNM-WNM interface, across which the 
density changes by a factor ~ 40 in ~ 1 pc. In addition to the 
large structures discussed above, there are small cold 'cloudlets' 
of size ~ 1 pc strung along field lines near the bottom of the 
panel. These resemble the cloudlets discovered by iHeilej il967t) 
and further discussed in iHeiles & Tro land d2003h . Figure^ shows 
a typical cloudlet; the local magnetic field lines are also shown. The 
maximum number density of this cloudlet is ~ 10 cm -3 ; its major 
(minor) axis is ~ 2 (1.2) pc, giving an axis ratio of ~ 5 : 3. The 
median column density is ~ 0.3 x 10 20 cm -2 . These numbers are 
in good quantitative agreement w ith the observations discussed in 
§ 8.3 of lHeiles & Troland fe003h . 



There is also a great deal of CNM away from the magnetic val- 
leys. In Figure(6] we show a cold Hi cloud elongated in the galactic 
plane at times (a) r = 7 and (b) r = 12. There are notable density 
variations and ripples along the CNM-WNM interface, where there 
is significant velocity shear. The major (minor) axis of the cloud at 
t = 12 is ~ 80 (1.35) pc, giving an axis ratio of ~ 60 : 1. The 
magnetic field makes an angle of < 22° with respect to the minor 
axis of the cloud. 

In light of the recent Hi surveys bv lHeiles & Troland Jiool) . 
it is of interest to calculate the fractions of interstellar mass in each 
thermal phase. In Figure [7] we show histograms of the fraction of 
the total mass in each thermal phase (CNM, unstable, and WNM) 
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Figure 7. Fraction of the total mass in each thermal phase (CNM, unstable, and WNM) at times r = — 13, —8,-3, 2, 7, and 12. Mass is binned by temperature 
in intervals of width A log[T/T re f (0)] = 0.1. The vertical dashed lines denote the temperature boundaries separating the different thermal phases. The total 
mass fraction in each of the phases is given as a percentage. The unit of time is C re f (0)/g = 6.05 Myr. 



at times r = — 13, —8, —3, 2, 7, and 12. Mass is binned by tem- 
perature in intervals of width A log[T/T re f (0)] = 0.1. The verti- 
cal dashed lines denote the^ temperature boundaries separating the 
different thermal phases (T m i n , T ma x)- The mass fraction in each 
phase is shown in each figure as a percentage. The reference state 
(r — —13) has all the mass in the WNM. As the Parker instabil- 
ity progresses, phase transitions are triggered, leading to a trans- 
fer of mass from the WNM, through the thermally-unstable phase, 



into the CNM. At the time when cold HI clouds have been formed 
and have settled into their final quasi-equilibrium states, the per- 
centages of mass in the CNM and WNM are approximately 40% 
and 60%, respectively. (This number is found by simply averag- 
ing the results between r = 7 and 12.) Enough time has elapsed 
for there to be an insignificant amount of mass in the thermally- 
unstable phase. 

Taken at face value, this calculation predicts the existence of 
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cold Hi clouds with maximum number density ~ 10 3 cm -3 , max- 
imum magnetic field strength ~ 20 fiG, and temperatures ~ 50 K. 
The mean values of the density and field strength in the cold clouds 
are found to be ~ 20 cm -3 and ~ 4.3 fiG, respectively. The rms 
value of the magnetic field strength in the CNM is ~ 6.4 fiG, while 
the maximum line-of-sight value of the field occurs for a line-of- 
sight in the azimuthal direction in the galactic plane and is ~ 13.5 
/j,G. The (two-dimensional) turbulent velocity dispersion inside the 
CNM clouds is ~ 1.6 km s~ ; this corresponds to a (thermal) Mach 
number ~ 2.8. The mass of a typical cold cloud is calculated by as- 
suming a thickness in the third (radial) direction corresponding to 
that of a galactic shock (~ 100 pc); we thus find the typical mass 
to be ~ 10 Mq . The spatial separation between cloud complexes 
is ~ 500 pc. An estimate of the maximum mass-to-fiux ratio of 
typical cold Hi clouds yields (M/$)/(M/<E>) crit =± 0.1 - 0.3, 
where (M/ $) cr it = (63G)~ 1//2 is the crit ical mass-to-flux ratio 
calculated bv lMouschovias & SpitzeJ dl976h . 



5 SUMMARY AND DISCUSSION 

We have performed two-dimensional simulations of the nonlinear 
evolution of the Parker instability properly accounting for phase 
transitions between the warm and cold neutral media. The evolu- 
tion is initiated by small-amplitude velocity perturbations that have, 
and at later times lead to, evolution that preserves the symmetry of 
the dominant midplane-crossing ('even') mode — the mode that 
allows field lines to cross the galactic plane. This mode has been 
shown by BMP97 to be naturally selected when starting from a 
spectrum of random velocity perturbations. The qualitative features 
of the nonlinear evolution are similar to those found by calculations 
that ignored phase transitions: arches of magnetic field lines rise 
high above the galactic plane, accompanied by accumulations of 
mass in the resulting magnetic valleys. However, there are impor- 
tant quantitative differences between the two sets of calculations. 
Phase transitions, a consequence of realistic heating and cooling 
included in the present calculation, lead to structures (cold clouds) 
with typical densities much greater and temperatures much smaller 
than those achieved by isothermal calculations, and which are in ex- 
cellent quantitative agreement with observations of cold Hi clouds. 

The main result of this work is that the Parker instability is 
directly responsible for the formation of interstellar clouds and 
cloud complexes, despite r ecent claims to the co ntrary by Kim et 
al. (1998, 2001, 2002) and lSantillan et all feOOCh . Cold Hi clouds 
are seen to form whose masses (~ 10 5 M©), mean magnetic field 
strengths (~ 4.3 fiG), rms magnetic field strengths (~ 6.4 fiG), 
mean densities (~ 20 cm -3 ), mass-to-fiux ratios (~ 0.1 — 0.3 
relative to critical), temperatures (~ 50 K), turbulent velocity dis- 
persions (~ 1.6 km s - ), and spatial separations (~ 500 pc) are all 
in agreement with observations. The maximum number density is 
~ 10 3 cm -3 and the maximum magnetic field strength ~ 20 fiG. 
These physical conditions are such that, had we included the addi- 
tional phase transition from atomic to molecular hydrogen, molec- 
ular clouds would have formed very rapidly within these giant Hi 
clouds. 

We find that the fraction of ISM mass in the WNM is ap- 
proximately 60% , in ex cellent agreement with observations by 
IHeiles & Trolandl J2003I) . The bulk of the CNM has a sheet-like 
morphology. We also note the presence of cold 'cloudlets' of size 
~ 1 pc strung along field lines. These cloudlets are both qual- 
itatively and quantitatively similar to those first discovered by 



Heilcs (1967). The median column density of a typical cloudlet is 
~ 0.3 x 10 20 cm -2 . 

An appare nt discrepancy between th is work and the observa- 
tional results of IHeiles & Trolana ( 120030 is the amount of WNM 
in the thermally-unstable phase. In the simulation, the maximum 
percentage of WNM in the unstable pha se is ~ 10%. By contras t, 
the observed percentage is ~ 48% by IHeiles & Trolan J d2003h . 
The discrepancy is easily understood as a consequence of the fact 
that the observations see the ISM as it is today, after clouds not 
only have formed but also some of them have given birth to stars, 
including OB associations, which alter the physical conditions of 
their parent clouds as well as the amount of matter in the warm un- 
stable phase. In other words, to account for the observed amount 
of gas in the thermally unstable phase, a calculation must include 
feedback from at least some of the energetic events that occur on 
relatively short timescales after clouds and stars have formed (e.g., 
stellar winds, UV radiation, or supernovae). Including such effects 
is beyond the scope of the present study, whose focus is cloud for- 
mation, not events that follow star formation. 

Our results differ subst antially from those of the r ecent three- 
dimensional simulations bv lKosinski & Hanaszl ( [2007h . Important 
differences also exist in the formulation of the problem. Their ne- 
glect of thermal conductivity, together with their very low resolu- 
tion (60 x 260 x 100 compared to our 2048 x 4096), results in 
thermally-unstable modes showing unresolved growth on their grid 
size of 10 pc (our grid resolution is ~ 0.2 pc). In addition, their as- 
sumed symmetry about the galactic plane limits the problem to the 
'odd' mode, which is not favored during the nonlinear evolution 
over the 'even' mode we have followed. Their heating and cool- 
ing functions imply a stable CNM phase at an unrealistically-high 
temperature of ~ 3000 — 4000 K (depending on the model). For 
the realistic heating and cooling functions we have employed, these 
temperatures lie in the thermally-unstable regime; we find that the 
thermally stable CNM phase exists below a temperature of 185 K. 

The results of this work lend strong support and quantify the 
scenario of cloud formation envisioned early on by Mouschovias 
et al. (1974). Taken in conjunction with the many successes of 
the theory of magnetic support of self-gravitating clouds and the 
ambipolar-diffusion-initiated fragmentation and star formation in 
such clouds (e.g., see § 4.2 of Mouschovias et al. 2006), the present 
work shows that magnetic fields play a crucial role not only in the 
formation of stars, but also in the formation and early evolution of 
their birthplaces, the interstellar clouds. 
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